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THEORETICAL STUDY OF REACTIVE AND NONRE ACTIVE 
TURBULENT COAXIAL JETS 

Roop N. Gupta* and N. T. Wakelyn 
Langley Research Center 

SUMMARY 

The hydrodynamic properties and the reaction kinetics of axisymmetric coaxial 
turbulent jets having steady mean quantities are investigated in this work. From the 
analysis, limited to free turbulent boundary -layer mixing of such jets, it is found that the 
two -equation model of turbulence is adequate for most nonreactive flows. For the reac- 
tive flows, where an allowance must be made for second -order correlations of concen- 
tration fluctuations in the finite -rate chemistry for initially inhomogeneous mixture, an 
equation similar to the concentration fluctuation equation of Spalding's model is suggested. 
This equation is much simpler to use than a similar proposal by Donaldson. For diffusion 
limited reactions, the eddy -breakup model based on concentration fluctuations is found 
satisfactory and simple to use. The theoretical results obtained from these various 
models are compared with some of the available experimental data. 

INTRODUCTION 

There are many practical applications of free turbulent mixing flow phenomena. 

For instance, the design of combustors for advanced aircraft (where weight and length 
should be minimized without the loss of combustion efficiency), jet ejectors, wakes behind 
airplanes and submarines, cooling water and waste dispersal in rivers, chimney plumes, 
and jet streams in the atmosphere are a few examples where the role played by free tur- 
bulent mixing is very important. Attention in this study is restricted to free turbulent 
boundary -layer mixing of reactive and nonreactive coaxial jets, even though many of the 
"real world" examples do not fit this description. The word "free" is used for the flows 
remote from walls. The word "boundary layer" is applied to the flow regions where a 
single predominant direction of flow exists and in which shear stresses, heat fluxes, and 
diffusion fluxes are significant only in directions at right angles to the predominant 
direction. 


*NRC-NASA Resident Research Associate, now at Aeronautical Engineering 
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Because of the practical importance of free turbulent boundary layers, a method for 
reliable prediction of these flows is greatly needed. Since Reynolds (ref. 1) formulated 
the equations of turbulent flow almost a century ago, various investigators have attempted 
to construct a model of shear flows based on a closure of these equations. The general 
approach, involving a modeling of appropriate terms in the equations for second -order 
correlations of fluctuating quantities, has been pursued in detail by Rotta (ref. 2) among 
others. More recently Donaldson and his colleagues (ref. 3) have attempted to develop a 
second-order closure model that can handle with one invariant model the classical as 
well as nonclassical shear layer problems (examples of the latter are generation of tur- 
bulence by the atmosphere and dispersal of chemically active species in the Earth' s 
atmosphere). They attempted in references 4 and 5 to couple their invariant diffusion 
model to the chemistry model. In the course of the analysis, they point out the situations 
for which second -order correlations involving concentration fluctuations are important 
and must be accounted for. In cases where the reaction rates are fast and the scale of 
turbulence is large, the reaction model based on mean-value chemistry can be substan- 
tially in error. An extreme case would be the flow of material by a given point such that 
large scale discrete parcels of the reactants alternately pass the point. For such a case, 
half of the time the flow is all one reactant and the other half of the time it is all the other 
reactant. In this admittedly hypothetical case, little or no reaction is possible, although 
mean-value chemistry will predict one. For such a poorly mixed system, the rate of 
removal of the reactants from the flow will be only weakly influenced by reaction rates 
and will instead be limited by molecular diffusion. On the other hand, if the scale of tur- 
bulence is small enough and the reaction rates are slow enough, diffusion is always fast 
enough to keep the reactants well mixed. For this case the depletion of reactants will, 
of course, be governed by reaction rates. In many laboratory flows, the diffusive scale 
of turbulence is small; hence, the mean-value chemistry model may be employed. How- 
ever, if the laboratory experiment is increased in size and the rest of the parameters 
held constant, the character of the flow will change because the scale of turbulence 
increases linearly with the size of the apparatus (ref. 6). 

Therefore, it is clear that for chemically active turbulent flows one has to be very 
careful in treating the second-order correlations involving the concentration fluctuation 
^cjcj^. However, in searching for a coupled turbulent diffusion and chemistry model 
which also allows for treatment of the concentration fluctuation or the "mixedness" 
term c^Cj, one has to look for a model which is acceptable in terms of the required 
accuracy as well as the effort and expense involved. The model proposed by Donaldson 
and his colleagues does not yet meet these requirements. Moreover, as they have also 
pointed out, the results obtained with their model should be used with caution because the 
length scale A employed in the calculation has a constant value. Their published results 
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so far are only pertinent where the equilibrium turbulence field is coupled to the chemis- 
try model. In view of these limitations and because of the complexity of their closure 
technique, an alternate approach with reasonable simplicity of computation and, as will be 
shown, adequate physics of the problem has been developed in this paper. 

In his study on turbulent flows with chemical reactions, Borghi (ref. 7) clearly 
points out the need for the modification of the molecular reaction rates by the concentra- 
tion fluctuations. However, his approach also has the limitation similar to that of refer- 
ences 3 to 5, namely, the complexity of closure technique for second -order moments. 

Also, a simple mixing-length hypothesis for the turbulence model has been employed in 
the analysis. 

Since the model presented herein is based on Spalding's work (ref. 8), his contri- 
bution will be described in detail. The Spalding three -equation (E k Wg) model, which 
solves the transport equations for turbulent kinetic energy E k , variance of the fluctua- 
tions of the vorticity W of the fluid (which may also be regarded as the length scale 
equation), and the variance of the fluctuating component of the mass fraction of the 
injected reactant g, appears to meet the requirements outlined earlier with adequate 
representation of the turbulent flow field. These equations are solved simultaneously 
with those governing the mean flow behavior. At a conference devoted to comparison of 
the predictions of various models (for nonreactive free turbulent shear flows) with each 
other and with experiment, Launder, Morse, Rodi, and Spalding (ref. 9) presented the 
results of the energy-dissipation ^E k e) model. Here e represents dissipation rate of 
the turbulence kinetic energy. They conclude that the E k e model is surpassed only by 
more complex "Reynolds stress" models, which are still not completely developed. In 
the E k e model the two transport equations to be solved are for E k and e = 
whereas in the E k W model the two turbulence quantities employed are E^ and 
W = E k /z^. (The symbol l indicates the length scale of turbulence.) There is little 
difference in the results using these two models for boundary -layer flows. For general 
flows, however, the E k e model is preferred (ref. 10). 

For round turbulent free jets, Spalding (ref. 8) has added a third equation for con- 
centration fluctuations to his E k W model. With this three -equation E k Wg model, 
Spalding has predicted concentration fluctuations that are in good agreement with experi- 
mental data for steady injection into a reservoir containing stagnant fluid of equal density. 

In the present treatment, a fourth equation, analogous to Spalding's equation for the 
variance of fluctuating component of the mass fraction of injected fluid, is proposed and 
used to predict the mixedness term cjcL It is thought that the present four -equation 
i.e., the Spalding three -equation model plus the equation for the mixedness 
term cjcj j may have the potential to handle adequately reacting turbulent flows for most 
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cases of interest. Also, the equations involved are simpler and fewer as compared with 
the other existing techniques which take into account the mixedness term c|cj (the 
approach of ref. 4, for example). Some predictions made by using the present model are 
given and compared with experimental data. 

For comparison with the proposed four -equation model, some results are also 
given by employing Prandtl's hypothesis (ref. 11) for free turbulent flows with mean- 
value chemistry. For fast chemical reactions and situations where the mixture is mainly 
composed of fragments of reactants and products, Mason and Spalding (ref. 12) have also 
suggested a turbulence controlled reaction model known as the eddy -breakup (EBU) model 
where the eddy breakup overrides chemical influence. The results from this model are 
included for nonpremixed reactants. 

With the inclusion of the length scale equation, the four -equation model is capable 
of handling reacting flows with nonequilibrium turbulence. The invariant coupled diffu- 
sion and chemistry model (ref. 4) is applicable only to equilibrium turbulence fields. 

The principal objectives of this investigation are briefly summarized as follows: 

(1) A four -equation model for free turbulence boundary -layer mixing of reactive 
coaxial jets is developed and tested. Since the second-order correlations involving con- 
centration fluctuations are important and must be accounted for in the calculation of cer- 
tain reactive flows, an equation for c!cj has been included in the model along with the 

transport equations for E^, W, and c'^. 

(2) The existing eddy -breakup reaction model employed with the Ej^W turbulence 
model is used to test the present four -equation model in the fast reaction limit. 

(3) Three versions of Prandtl's hypothesis for free turbulent flows are employed 
with mean-value chemistry for comparison with the models given in (1) and (2). 

The present analysis has been limited to one-way forward reaction to keep the 
model and the calculation procedure simple. For this purpose the following reaction 
with mild heat release is considered: 

N0 2 C1 + NO — * NOC1 + N0 2 

Since this reaction goes to at least 99 percent completion at about room tempera- 
ture, the fluctuations in temperature or the reaction rate coefficient are not considered 
to be significant. The extension of the present work for a more complex reaction system 
with large heat release is likely to require considerable effort. Nevertheless, efforts 
similar to the present work should prove more practical and appealing, especially if 
adequate physics of turbulence is to be accounted for. 
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SYMBOLS 


a jet nozzle radius, meters 

C constant used in various equations 

c mass fraction of species 

Cp specific heat 

Djg binary coefficient of diffusion 

d jet nozzle diameter, meters 

Ejj turbulent kinetic energy 

g mean -square concentration fluctuation 

H stagnation enthalpy 

h =2 C i h i 

i 

hq enthalpy of ith species 

J turbulent flux tensor 

K thermal conductivity of mixture 

k mean reaction rate coefficient 

k' fluctuating component of reaction rate coefficient 

l local average length scale of turbulence, see equation (12) 

M molecular weight of species 

M m ix molecular weight of mixture 

m” rate of mass transfer across outer boundary 



N L e 

molecular Lewis number, 

N Le,t 

turbulent Lewis number, ^ 

Npr 

molecular Prandtl number, 

N Pr,t 

turbulent Prandtl number, 

N Re 

Reynolds number 

N Sc 

molecular Schmidt number, 

N Sc,t 

turbulent Schmidt number, 


p 

R 


R i 


ci 


S 

St 


>U 


T 

t 

U 


K 


CpM 

K 


9Cj 

3y/ ay 


N 


Pr 


N Le 

N Pr,t 

N L e,t 

pressure 

universal gas constant 

rate of production of ith species 

radial coordinate 


Cj <t + Cj p 

concentration half -radius where Cj = ? ? 


velocity half -radius where u = 


u t + u e 


generation-dissipation term 

downstream component of generation -dissipation term 
upstream component of generation-dissipation term 
temperature 
time 

mean velocity 
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I 


u mean axial velocity 

u' fluctuating component of axial velocity 

v mean radial velocity 

v' fluctuating component of radial velocity 

W turbulent vorticity fluctuation 

w' fluctuating component of tangential velocity 

x axial coordinate 

y vertical coordinate 

a constant for external control of eddy viscosity, equations (31) to (36) 

6x forward step size in numerical solution 

k eddy viscosity coefficient, equations (32) and (34) 

i< eddy viscosity coefficient for Donaldson -Gray model, equations (35) and (36) 

A length scale, equations (23) to (29) 

X microscale length, equations (27) to (29) 

ju laminar viscosity 

/it turbulent viscosity, -uVy/|^ 

v kinematic viscosity 

p density 

ct parameter in finite -difference equation 

r shear stress 
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0 


general dependent variable 


X mole fraction 

4/ stream function 

co dimensionless stream function 

Subscripts: 

A species A 

B species B 

chem chemistry 

t center line 

Ek relating to turbulent kinetic energy 

e external stream 

ebu relating to eddy -breakup model 

eff effective 

i species i 

j species j 

max maximum 

o nozzle exit plane 

t turbulent flow quantity 

W relating to turbulent vorticity fluctuation 

(p relating to general dependent variable cp 
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II 


Superscripts: 

time mean 
’ fluctuating 

Chemical symbols: 

NO nitric oxide 

NOC1 nitrosyl chloride 

NO 2 nitrogen dioxide 

NO 2 CI nitryl chloride 


ANALYSIS 


Time -Averaged Boundary -Layer Equations for Reactive 
Axisymmetric Turbulent Flow 

The following conservation equations (refs, 13, 14, and 15) govern free shear layer 
turbulent mixing of coflowing axisymmetric streams (fig. 1) undergoing chemical 
reactions: 

Momentum: 


9u 3u 9p 1 8 
ou — + ov — = — — + — — 
H dx H dr dx r dr 


(<* + "t)r # 


U) 


Species: 

acj 


9c i 1 a 


Oil T — + PV „ 

^ dx ^ dr r dr 


Energy: 


pu |H +pv |H = IJ. 

dx dr r dr 


M t \ 

+ r 


N Sc N Sc,t/ 


iL_ +J±J\ r 9 H 


+ R i 


(2) 


N Pr N Pr,t/" 9r 


+ 1 f 

r dr 


M ( X ' N Pr) + ' Npr.t) 


_a_ Hi 
dr 2 


r ar 


^Le ~ * ^Le,t “ ^ V u 9c i 

M N Pr + Mt Np r t j r Z, h i ar 


(3) 
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Mass: 


9(pu) , 1 9(pvr) _ 0 
9x r ar 


(4) 


In addition, for ideal reactive gases the following equation of state is also used: 


P = P 



(5) 


where R is the universal gas constant and M mix is the molecular weight of mixture. 
The density fluctuations of individual species have been neglected. 

The production term in equation (2) may be expressed in terms of time- 
average and fluctuating quantities: 



where it is assumed that k[ = 0. This assumption is true for an isothermal system. 

The following initial and boundary conditions are used to solve the conservation 
equations (1) to (4): 


X = 0: 

u = u(r) 

c { = C i (r) 

H = H(r) 

r = 0: 

au _ aH _ 

ac4 

— = 0 



ar ~ ar 

ar 


r oo: 

u - U e 

H-H e 

c i — c i,e 


For axisymmetric jets, the value of the turbulent Prandtl and Schmidt numbers is 
approximately 0.7 (ref. 15); this implies that the turbulent Lewis number may be taken 
as unity. By assuming the laminar Lewis number to be unity also, the last term in the 
energy equation (3) disappears. Also at low subsonic speeds, the second term on the 
right-hand side of equation (3) containing velocity is negligible. Thus, energy equa- 
tion (3) may be simplified to 


pu|» + pv|5 
ax ar 


I JL 

/ M M t ) 

r 8H 

r ar 

\Np r N Pr>t ) 

r ar 


(8) 
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Further, in order to keep the chemistry simple, an irreversible one -step reaction 
of the type 


A + B 


C + D 


0) 


has been used in the computations. 

For obtaining the laminar viscosity p of gas mixtures, the following semi- 
empirical formula of Wilke (ref. 16) is used: 

(10) 


(ID 

To complete the description of the problem, the eddy viscosity and concentration 
fluctuations must be specified. In the following sections several models are discussed. 

Two -Equation (E k W) Model of Turbulence 

The turbulence model described was developed by Spalding and his colleagues 
(refs. 8, 9, 10, 17, 18, and 19). It employs two parabolic partial differential equations 
as follows: 

(1) One equation governing the distribution through the turbulent field of the specific 

kinetic energy of the fluctuating motion 

(2) A second equation governing the distribution of a property W having the 

dimension of frequency squared 

E k 

The quantity W = — may be interpreted as the time -mean square (or variance) 
l 2 

of the vorticity fluctuations. In fact, the equation employed can be for l (local^average 
length scale of the turbulence), E k Z, or in place of the equation for = W. 

The turbulent effective viscosity is related to E k and l by 

M t = pE k 1/2 J (12) 


where 




Xi^ 


i 

j 


1 ( Mi 1/2 

*«* 7j( 1+ 5j 


1 + 


'PA1/2/MA1/4 


Mi, 
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with 


E k = 



+ w' 2 ) 


(13) 


l = 



(14) 


The equations governing E k and W for low subsonic flows are, respectively, 


pu 


9E k 9E k 

i 


Advection 


= ii-lr 


H 


aE 


r9r \ N Pr,Ek 9r 
Diffusion 


r) + P E k 1/2 *(f7) 2 ' C D p 


E k 


3/2 


Generation 


Dissipation 


(15) 


aw aw 
pu !T + pv lF 


Advection 


- -r-|r 


H 


r 9r \ N Sc,W 9r 

Diffusion 


9 W\ p p 1 /2 7 ( d U j p 

— + C lP E k Z — s + C 3 p 


. 2„\ 2 _ EJ / 2 


.ar* 


Generation 


_ E.- 1 / 2 


- C 2P' 


-W 


Decay 


(16) 


In these equations Np r g k and Ng c are Prandtl and Schmidt numbers, respectively. 
Further, equations (15) and (16) are solved subject to the following initial and boundary 
conditions: 


x = 0: 


E k = E k( r ) 


W = W(r) 


r = 0: 


££k = aw = 0 

ar ar 


(17) 


r - E k - E k,e 


W-W e J 


The values of the various constants appearing in equations (15) and (16) are taken from 
reference 16 and are given as follows: 
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Ill 


C D = 0.09 
C x = 3.50 
C 2 = 0.17 
C 3 = 1.48 

N Pr, Ek = °- 90 
N Sc,W = °- 90 








(18) 


Equation (16) has been derived (ref. 20) by analogy to equation (15). The turbulent energy 
generation term is » the shear stress, r = The eddy viscosity is taken 

equal to pE^^l. The dissipation is assumed proportional to The constant 

Cj) is assigned a value of 0.09 to satisfy the energy equilibrium characteristics of 
constant stress wall layers. 

Though the equations employed in the turbulence model are for E^ and W |same 
as Ej other models solve for different combinations of E^ and l. For example, 
equations for e = Ej or E^Z have been employed in place of the equation for W 
(refs. 9 and 10). Whereas there is little difference in the results between these two 
models for boundary -layer flows, recent results (ref. 3) indicate that for general flows 
the e equation is preferable to the W and E^ equations. 

Correlations of Concentration Fluctuations 

A third equation in the turbulence model developed by Spalding (ref. 8) to predict 
concentration fluctuations is given as 


9gAA 
pu — — 


+ pv 


8 SAA 

3r 


1 J_L ag AA\ 

r 9r \ N Sc,g 9r 1 


+ 


C g lP E k 



E k 1/2 

< -'g2 P l &AA 


(19) 


where 

_ <2 ~ ~ 
g AA = C A = C A C A 
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Here Ng c g has the significance of "effective Schmidt number" and has been assigned a 
value of 0.7 (ref. 8). The constants Cgj and Cg 2 are given values of 2.7 and 0.134, 
respectively, from reference 8. 

A fourth equation similar to equation (19) is proposed for §ab( c A c b) : 


pu 


3gAB 

9x 


+ pv 


SgAB 

ar 


Advection 


la/ Mt 9gAB\ 
r 9r \ N Sc,g dr ) 

Diffusion 



Generation 


Ek 1/2 

c g 2 f — — sab 

Decay 

(20) 


The following initial and boundary conditions have been employed to solve equa- 
tions (19) and (20): 


x = 0: 
r = 0: 


r 


oo; 


gAA gAA^ r ^ 

f£AA = f£AB 
ar ar 

&AA (Saa) 6 


Sab - SAB^ r ^ 

> 

sab • (sab ) 6 


( 21 ) 


The values used for Cgj and Cg£ appearing in equation (20) are the same as 
those suggested for equation (19). These values might need some change for certain 
reactions, since Cg^ and Cg 2 have been chosen in reference 8 to obtain agreement 
with the experimental data on fluctuations of concentration of the injected fluid. However, 
it may be noted that Cgj has a value close to that of Cj whereas Cg 2 is closer 
to C£. This is to be expected in light of the similar significances of the relevant terms 
in the equations. 

In order to compare equations (19) and (20) with those given in reference 21 
(pp. 7-1 to 7-20), these equations are further simplified for the case of equilibrium 
turbulence. Thus, by assuming that generation of g equals dissipation of g locally, 
the following equations are obtained from equations (19) and (20): 


_ = 2.7 ; 2/££A' 

Saa 0.134 l ar 


C A C A = 20. 1491“ 


2 dc A 


ar 


(22a) 
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I 


2 ( dc A ( dc B 


c^- 20. 1491' ^ -^ 


(22b) 


Equations (22a) and (22b) compare directly with the following equations given by 
Donaldson (ref. 21) for superequilibrium theory: 


C A C A = CCa2 ( w 


(23) 




(24) 


where 


CC = 1.706 

Equations (22a) and (22b) are equivalent to equations (23) and (24) if the length scale Z 
used by Spalding is related to the length scale A used by Donaldson by the relation 

Z = 0.291 A (25) 


or, in terms of Spalding's coefficient, 


Z = c d 1/2 a 


(26) 


The chemical mixedness term cjcj appearing in equation (6) is evaluated by using 
equation (20) rather than its simplified form (eq. (22b)).- 

Equations (19) and (20) can be further compared with those obtained in reference 4 
for the invariant coupled diffusion and chemistry model. The authors of reference 4 have 
obtained equations of the following form for g^g and g^^ * or an atmospheric shear 
layer with constant density: 


u 


9g AB 

ax 


_ — — V 9c A 9 

“ " C A V ~9y" " C B V "9y 


T — t ac B 


+ i|( a *nK 


9 §Ab'\ 9 g AB 

+ V c\ 


ay 


ay £ 


g AD 

- 2v —2 ' k Al c A g BB + c B g AB + c 


, , 2 \ 
A C B / 


" k B\ c A g AB + c B g AA + C A 


,2 , ] 

A C B] 


(27) 
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u !!aa = + i( A *^ k !£aa) + , % 

A 9 y 3y\ v K ay / ay 2 


ax 


g AA 


2V “72" " 2k A\ c A g AB + c B g AA + C A C B 


(28) 


where y is the vertical distance, v' is the fluctuating component of the vertical 
velocity, A* is the scale associated with velocity diffusion terms, and X is the dissi- 
pative scale. The quantity A which appears in the superequilibrium approximation 
(eqs. (23) and (24)) is the isotropy scale. In general, there will also be a scale A** 
associated with the pressure diffusion terms in the invariant turbulence model. In a 
simple model all three A's may be taken as equal; the following relationship may be 
used to relate them to X (refs. 21 and 22): 


X = ^- v - - (29) 

V a +bN Re,A 

where a and b are constants and Nj^ e ^ is a Reynolds number based on A and E^. 

In the analysis presented in reference 4, the third-order correlation of equations (27) 
and (28) has been modeled in terms of the second-order correlations. However, the 
results are applicable only to the case where the chemistry model is coupled to the equi- 
librium field of turbulence. The length scale A is prescribed as an algebraic value. 

Also, extra equations are needed for modeling c^v' and CgV* and further assumptions 
are made to model some of the higher -order terms appearing in these equations. In view 
of these assumptions and the complexity of their turbulence model, it is felt that equa- 
tions (19) and (20) may be employed with lesser computational effort. Moreover, equa- 
tion (19) has been used to predict (ref. 8) concentration fluctuations with good agreement 
with the experimental data for the case of steady injection into a reservoir containing 
stagnant fluid of equal density. 

Eddy -Breakup (EBU) Model for Reaction Rates 

For fast chemical reactions, a simpler model which adequately accounts for the 
fluctuations in temperature and concentrations has been proposed by Spalding (ref. 12). 

The model known as the eddy -breakup model involves the calculation of g^^ either 
from its differential equation (eq. (19)) or its approximate algebraic form (eq. (22a)). In 
this model the reactions are essentially "diffusion limited;" that is, the reaction rate is 
determined by the rate at which one of the reactants is supplied by diffusion rather than 
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by chemical kinetics. Spalding has proposed the following eddy -breakup reaction rate 
expression for premixed turbulent flames of high temperature: 


1/2 E k 1//2 

^i ^ebu^AA i 


(30) 


where C e ^ u is a universal constant for high Reynolds number flow and has an approxi- 
mate value of 0.53. Equation (30) is employed to model the reaction between NO 2 CI 
and NO, which is a quite fast reaction with mild heat release (ref. 23). For the case 
analyzed here, NO 2 CI and NO are not premixed. The rate prediction depends upon the 
concentration fluctuations of the injected reactant (n 0 2 C1). However, to satisfy the 
requirement that no reaction is predicted whenever one of the two reactants is zero, the 
smaller of the NO 2 CI or NO concentration fluctuations is employed in the EBU model. 
Comparison of the results is made with those for the case where the reaction rates are 
predicted from expressions (6) and (20). The turbulence model employed with these two 
chemistry models is the Spalding E k W model of turbulence. 


Turbulence Models Employing Prandtl's Hypothesis 

The following eddy viscosity models based on Prandtl’s hypothesis for free turbulent 
flows were also used to predict the axisymmetric turbulent jet velocity profiles. 

(1) The Ting-Libby model (ref. 24): 

Initial region: 


Upstream of axial position where 
li t = pe = a-(0.00137)|u o 


Uf u e 

u 0 - u e 




Developed region: 

M t = pe = oud 1/2 |U t - U e |p^j ($f 

where 



(31) 


(32) 
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' 1/2 


value of 7 ] where u = 


U<t + u e 


eddy diffusivity 


k eddy viscosity coefficient given by Schlichting (ref. 10) as 0.025 

Here a may be used to change the value of k. 

(2) The Ferri model (ref. 25): 

Initial region: 

M t = pe= a(0.00137)| P<t U t -p e U e | 

Developed region: 

p t = pe= Q!/cb 1/2 |p t U t -p e U e | 

In equation (34), bjy 2 is the value of r where 
pu p t Ut + P e U e 


(33) 


(34) 


(3) The Donaldson-Gray model (ref. 26): 
Initial region: 


H = P e = ai <( r l/2 ~ r i )p 


U f U e 


(35) 


For subsonic flows, k = 0.0468. In equation (35), r j / 2 is the value of r where 

u <t - u e 

u = and r^ is the inner mixing-zone radius defined as the value of r where 


u - U e 

U 0 - U e 


= 0.95 


Developed region: 


u f u e 


li t = pe = a K v 1/2 p ^ 


(36) 
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These three eddy viscosity models predict the potential core and the axial decay 
with very little difference. In this work, however, the Donaldson -Gray model is used for 
comparisons with the two -equation model. All these models are, of course, coupled to 
the chemistry models in which reaction rates are specified by simple time averaging of 
the instantaneous Arrhenius expression with neglect of the mixedness term cjcj and 
fluctuations in the reaction rate coefficient (k'). 

In summary, the Spalding E k Wg model (ref. 8 ) is investigated in detail. An extra 
equation analogous to Spalding's concentration fluctuation equation is added to this model 
to evaluate the mixedness term c{cj. This along with the mean-value term is then used 
to evaluate the reaction rate expression from the Arrhenius rate equation. The eddy- 
breakup reaction model with coupling to the E k Wg model and the eddy -viscosity model 
based on Prandtl' s hypothesis with mean -value chemistry are included for comparison 
only. The comparison would reflect the magnitude of error involved when Prandtl's 
hypothesis is employed for certain predictions. 


Transformation of Equations 

For the axisymmetric jet problem considered herein, it is convenient to use 
Von Mises transformation. The stream function can be defined from the continuity 
equation (eq. (4)) by its partial derivatives: 


pur = 

K dr 


pvr 


dip 

dx 




(37) 


All the differential equations of interest are parabolic partial differential equations of the 
form 


pu !£ + pv !? “ r £{ rJ <p) + S <p 


(38) 


where 


T _ Meff 9(f> 

J 4> = T^ dr 

By applying transformation (37) to equation (38), the following equation is obtained: 


M = »L i^ff r 2 acA + J_ g 
ax ai pr <j 0 ai jyj pu <P 


(39) 
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When 




pur dr 


pur dr 


is introduced and 


9^e 

ax 


" r e m 


ft 

e 


from the definition of \p, then equation (39) may be written as 

9(f) r e I ^ l e 90 _ 9 /r^pu ^eff 90\ , J_ g 
9x + \p e w dco ~ 9w^ e 2 9c oj pu 0 


(40) 


(41) 


(42) 


where rhg is the rate of mass transfer across the e boundary. (See fig. 2.) The 
introduction of the new variable o>, defined by expression (40), helps to keep the flow 
field between 0 and 1 and the number of computational points well distributed. In equa- 
tion (42), the terms on the left-hand side are representative of advection by the mean flow 
and those on the right-hand side represent the diffusion and source/sink terms, respec- 
tively, of the entity 0. 

Equation (42), with the reaction models given by equations (6) and (30), is solved 
with the help of the Patankar -Spalding GENMIX finite -difference procedure (ref. 13) with 
some modifications suggested in reference 27. Some of the changes are described 
briefly subsequently. 

Predictions with the turbulent eddy -viscosity models based on Prandtl's hypothesis 
and the mean-value chemistry are made by employing the method of calculation of refer- 
ence 28. It is a mixed implicit /explicit finite -difference scheme. In this scheme, the 
stability problems inherent in fully explicit finite -difference schemes are shown to be 
eliminated, and stable step sizes are shown to be increased by up to 4 orders of magni- 
tude. Details of the method of calculation can be found in reference 28. 


Selection of Chemicals 

As pointed out in the beginning, a one -step reaction of the type A + B — ► C + D 
has been used in the analysis to keep the chemistry simple. A reaction which falls in 
this category is the one between NO 2 CI and NO. The stoichiometric reaction between 
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NC^Cl and NO, which goes to 99 percent completion at room temperature, is believed to 
proceed by way of the following one-step reaction (refs. 23, 29, and 30): 


k 

N0 2 C1 + NO — - N0 2 + NOC1 


(43) 


where k = 1.39552 x lO - -^ cm^/molecule-sec at 300 K (ref. 23 *). Reference 29 points 
out that the forward reaction in this case is faster than the backward reaction by about 
4 orders of magnitude. The progress of the reaction in equation (43) can therefore be 
followed by calculation or measurement of the concentration of a single species. 

The four -equation theoretical model proposed in this work to allow for the mixed - 
ness term c[cj, however, is not limited to irreversible bimolecular reactions. Reac- 
tion (43) has been chosen with a view to keep both the theoretical predictions as well as 
the experimental measurements simple. Reversible and multistep reactions may be 
treated by extension of the four -equation model with appropriate equations, but the validity 
of such modeling would first need the verification of equation (20) from measurements. 

Initial Conditions 

The profiles of u, E k , and l needed to start the integration at the nozzle exit 
plane have not been measured in some of the experimental data included herein for com- 
parison. Therefore, for prediction of these flows, the initial conditions were assumed as 
follows. The velocity profiles were taken as uniform at the exit and the streams were 
assumed to carry no turbulence. The calculations were started from estimated distribu- 
tions of u, E k , and l in a very thin mixing layer at the edge of the nozzle lip. In the 
predictions, a self-similar mixing layer evolved. Therefore, the initial distributions of 
u, E k , and l in the very thin layer were found to have no appreciable influence on the 
downstream results. This approach was used by Rodi (ref. 31) in the prediction of non- 
reactive free turbulent boundary layers. 

Where initial profiles of velocity at the nozzle exit plane are available, the following 
relation between turbulent shear stress and turbulent kinetic energy is employed (ref. 32): 

T t = 0.3pE k (44) 


* The value of k given in reference 23 is incorrectly quoted in both references 29 
and 30. 
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or from equation (12), 


P E k 1/2 l |f = °* 3 P E k 


or 

F 1/2 _ Uq l s(u/U 0 ) 

^ 0.3 a 9(r/a) 


( 45 ) 


o 

To evaluate equation (45) at the nozzle exit plane, the length scale l is taken as 
the constant fraction C-^Crj of the shear layer width (ref. 33) and determined by using 
the following formula: 


1 = C 1 C 7|Yo. 1 " YO.ol (46) 

with 

U (yp.i) - u e _ 0 1 

U t - U e 

and 

U ( y °- 9 ) - U e = 0.9 

u t - U e 

In reference 33, values of 0.08 and 0.875 have been suggested for constants Cj 
and C 7, respectively. Typical initial velocity profiles and the characteristic shear 
widths for the inner and outer regions of flow are depicted in figure 2. 


Numerical Solution 

The Patankar- Spalding GENMIX finite -difference procedure (ref. 13) has been 
employed for the solution of equation (42), with the reaction models given by equations (6) 
and (30). Some modifications suggested in reference 27 have been adopted to improve 
the accuracy of this technique. These modifications are briefly described in this section. 
The calculations using the turbulent eddy viscosity models based on Prandtl’s hypothesis 
(eqs. (31) to (36)) and mean-value chemistry have been made by employing the program 
described in detail in reference 28. 

The essence of the GENMIX program (ref. 12) is the transformation of equations (1) 
to (3) and equations (15), (16), (19), and (20) to the independent variables (x,u>) with c o 


O 

In this analysis, the single length scale l has been assumed to influence the 
magnitudes of both the transport and the dissipation processes in a fully turbulent fluid. 
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defined by equation (40). The details of the numerical approximation are given in refer- 
ence 13, and reference 27 contains some of these details along with modifications. Some 
of these details are listed below: 

(1) Patankar and Spalding, using a marching integration procedure, formulated a 
finite -difference equation by integrating equation (42) around a microcontrol volume by 
using a trapezoidal method of integration. This linearization between grid points in the 
cross-stream direction can result in large truncation errors. Especially, when there is 
a large rate of increase in <p at one grid point, then flux conservation causes a decrease 
in cp at the neighboring grid point - a physically unreasonable situation. However, a 
better approximation of values between grid points can be made with Taylor series 
expressions to form the finite -difference equations than can be made by means of linear 
interpolation. 

(2) A variable parameter a is included in the difference equations so that a 
forward difference or a central difference scheme may be adopted. The Patankar - 
Spalding method employs an implicit forward difference method so that cr = 1. Formally, 
the truncation errors are minimum for the Crank-Nicholson central difference method 
for which a = 0.5. However, this may not always be true. For example, the finite - 
difference method for Tollmien's jet in still surroundings (ref. 34) was found to give best 
results (ref. 27) for large step sizes when a = 1. This is due to the fact that, for this 
problem, use of the upstream value of viscosity slightly overestimates the diffusion term, 
whereas the calculation of d^cp/d from the downstream stations results in a slight 
underestimation. These two effects tend to balance each other. Generally speaking, a 
value of a between 0.5 and 1 may be used to minimize the calculation time for a given 
problem. For the axisymmetric jet analyzed herein, cr = 0.5 was employed. 

(3) The decay/production term is represented in a way similar to that of the 
Patankar -Spalding method: 

d = pu S <£ = S U + aS D^ + (*• " ct ) s D < % ( 47 ) 

where the subscripts U and D represent upstream and downstream values, 
respectively. 

(4) As suggested by Patankar and Spalding, a forward step size equal to a small 
fraction of the mixing layer thickness is used. However, when the generation-dissipation 
term in equation (42) becomes large compared with the diffusion and advection terms, 
then the following approximate equation may be used: 


1 g 

ax " p5 


(48) 
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or, in finite -difference form, by using equation (47), 

0^1 + (1 - cr)Sp5xJ + SySx 
^ (1 - a)Sj)6x 


(49) 


Equation (49) yields the following limit for the forward step size: 


6x < 


1 

CT (Sl>) max 


(50) 


where ( s D) max is the l ar S est S D at a g iven cross section. For the present problem, 
(Sx)max ^ as b een taken as 



for the results to be accurate. For extremely fast reactions, some other measures 
should be taken so that the step size is not very small. 

Mention must now be made of the specification of entrainment rate to avoid large 
inaccuracies. The advantage of the Patankar- Spalding coordinate system is that the same 
number of grid points may be used to cover an increasing field width as the solution pro- 
ceeds downstream. This is accomplished by specifying the entrainment rate r e m" in 
equation (42) so as to increase the mass flow under consideration with each step. The 

^ .ft 

r e m e ad) 

entrainment rate determines the magnitude of the advection term — - — i>) — ^ in equa- 

tion (42) and, therefore, influences the whole flow field. In particular, advection should 
be minimal in the initial development of the solution if profiles are not smooth. At the 
same time another condition must also be met so that the finite -difference solution cor- 
rectly represents the physics of the problem; namely, at the outer free boundary of the 
flow, the gradient of <p must approach zero. These requirements on entrainment in the 
initial developing region of the jet can be met by providing a number of outer points with 
values close to the free-stream value. As the solution advances downstream, the jet 
spreads through these points and appreciable entrainment does not occur until significant 
variations appear at points adjacent to the outer boundary. 

The formulas used for specification of entrainment rate are arbitrary as long as 
they satisfy the previously mentioned conditions. The effect of the specified entrainment 
rate at each step is not known until the difference equation has been solved. Therefore, 
the process becomes one of guesswork to meet the earlier mentioned requirements. The 
following procedure has been adopted herein to specify the entrainment rate: It is first 
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assigned an arbitrary value. If the chosen value is too small (as exhibited by the appear- 
ance of a gradient at the outer edge), it is doubled. On the other hand, if the value used 
for the entrainment rate is too large (as exhibited by <p taking negative values), it is 
halved. 

The actual computation procedure employed in this work is outlined. After speci- 
fication of the initial profiles of the dependent variables and other parameters, the 
dimensionless cross -stream variable u> and all functions thereof are calculated. The 
numerical solution is now advanced step-by-step as follows: 

(a) Obtain radius r from cross -stream variable co and density 

(b) Calculate viscosity across stream from turbulence energy and length scale 

(c) Calculate decay/production terms of differential equations 

(d) Set forward step size and entrainment rate 

(e) Calculate finite -difference coefficients 

(f) Adjust boundary values and solve finite -difference equations by simple succes- 

sive substitution method 

RESULTS AND DISCUSSION 

The predictions made by employing the various turbulence models considered are 
discussed for reactive and nonreactive flows and are compared with some of the available 
experimental data. The predictions have been made for coaxial axisymmetric jets sub- 
ject to the initial conditions presented in table I. 


TABLE I.- INITIAL CONDITIONS FOR TEST CALCULATIONS 


Case 

Inner 

jet 

Outer 

jet 

u e/ u o 

U o> 

m/sec 

Nozzle 

diam., 

cm 

T e/ T o 

T e , 

K 

Inner jet 
density, 

kg/m ^ 

Outer jet 
density, 

kg/m^ 

1 

Krypton 

Air 

0.25 

36.580 

0.635 

1 

300 

3.4300 

1.18 

2 

Air 

Air 

.25 

36.580 

.635 

1 

300 

1.1800 

1.18 

3 

Hydrogen 

Air 

.33 

731.707 

1.270 

— 

— 

.0818 

1.18 

4 

no 2 ci 

NO 

.25 

36.580 

.635 

1 

300 

3.3300 

1.23 


Since krypton is an inert gas, the krypton -air mixing case (case 1) has been 
chosen to study the turbulent mixing in the absence of chemical reactions. Further, the 
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molecular weight of krypton approximately equals that of NO 2 CI and the molecular weight 
of air approximately equals that of NO. Therefore, a comparison between the nonreactive 
krypton -air mixing and the reactive NO 2 CI-NO mixing can tell us how the reactions can 
modify the mixing. 

Figures 3 and 4 show results for krypton-air mixing. These figures give the 
center -line velocity and composition decay for the turbulence models based on Prandtl's 
hypothesis and the Spalding two-equation model. The three turbulence models - namely, 
Ferri, Ting-Libby, and Donaldson-Gray - based on Prandtl's hypothesis predict the 
potential core and the axial decay with very little difference. Out of these three models, 
the Donaldson -Gray model is used for further comparisons. The Spalding two-equation 
model also predicts the potential core with almost no difference. However, this model 
predicts faster decay in the axial values of velocity and concentration. The Spalding 
model is preferred in this work because further comparisons with experimental data indi- 
cate this to be a better model. 

Figures 5 and 6 give the results for air-air mixing (case 2). The conditions for 
case 2 correspond to Forstall series E experiment (ref. 15). The initial boundary -layer 
profile of reference 35 is used. In figure 5(a) the experimental data of reference 15 are 
compared with the predictions made by using the Donaldson-Gray model. In this model a 
factor a is employed which may be used to change the value of the eddy viscosity coeffi- 
cient k. In addition to the usual a = 1.0, a power -law variation of a with x/d is 
also incorporated. Although some improvement in the prediction of center -line velocity 
decay is noticeable with the use of power-law variation of a, the results with a = 1 do 
not appear much different. Therefore, a will be used with a value of unity except where 
otherwise indicated. 

In figure 5(b) the results of the calculations from different turbulence models are 
shown. Clearly the predictions from Spalding's two -equation model are superior to the 
other two models when compared with the experimental data of Forstall and Shapiro 
(ref. 15). The predictions from the Donaldson-Hilst invariant model have been made by 
Hilst. The initial turbulence for the model has been specified in the same way as that 
used for the two -equation Spalding model, namely, 


K 1/2 ), 


initial 


h (0 - 084) 


8 (“/ u o)l 

9(r/a) 


initial 


The predictions with the invariant model appear unsatisfactory. It predicts almost no 
potential core and a very rapid decay as compared with the data of Forstall and Shapiro. 
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The prediction with the invariant model might be improved if some of the "constants" 
employed in the model are modified. 

Figure 6 gives some additional comparisons for the air -air mixing case. Here the 
various models are used to predict the velocity half -radius r u . Once again, the predic- 
tions from the two-equation model appear more satisfactory as compared with the invari- 
ant model or the Donaldson-Gray model. 

For the hydrogen-air mixing (case 3), the two-equation model predictions are 
excellent. (See figs. 7 and 8.) The conditions for case 3 are the ones used in refer- 
ence 36. In the initial mixing region ^3 1^ 15^, when the field of turbulence is non- 
equilibrium, the prediction methods (such as that of the Donaldson-Gray method) without 
a length scale equation appear to do a poor job. This result may be because, in eddy 
diffusivity theories of the type of Donaldson-Gray, the turbulent transport at a point is 
independent of the value of transport at another point. In the two -equation model, this 
deficiency seems to be removed by the use of a length scale equation. For this reason 
the Spalding two -equation model is considered better. 

The predictions for the reactive mixing of NO 2 CI-NO coaxial jets (case 4) are con- 
tained in figures 9 to 11. Figure 9 shows the axial decay of the center -line velocity and 
the longitudinal turbulence intensity. The Donaldson -Gray and the two-equation model of 
turbulence both predict the center -line velocity decay almost in the same way. Also, the 
inclusion or omission of the mixedness term does not influence the prediction of velocity 
decay. However, the mixedness term does play a significant role in the prediction of 
center -line composition decay as shown in figure 10. The Donaldson-Gray model, which 
employs the mean-value chemistry, predicts complete depletion of NO 2 CI by ^ ~ 14. 

The two -equation model of turbulence, with the inclusion of mixedness term in the chem- 
istry model, predicts that NO 2 CI still has a value of about 20 percent of its nozzle concen- 
tration at £ « 14. Therefore, a reaction is predicted beyond this axial location. The 
eddy -breakup model, which employs the expression for the eddy -breakup reaction rate to 
predict the diffusion limited flows with concentration fluctuations, predicts results com- 
parable to the Spalding two -equation model employing the chemistry model with mixed- 
ness term. The mixedness terms are predicted from equation (20) which is quite simple 
to employ with the two -equation model of turbulence. Some accurate and careful mea- 
surements are, however, desirable to firmly establish equation (20). 

The importance of retaining the mixedness term is more evident in figure 11. 
Without the inclusion of this term, the depletion rate of NO 2 CI is overpredicted by about 
100 percent. Therefore, this term should be included in a chemistry model and compari- 
son with experimental data is required to assess its importance. 
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CONCLUDING REMARKS 


This investigation is concerned with the hydrodynamic properties and the reaction 
kinetics of axisymmetric coaxial turbulent jets. The analysis is limited to the free turbu- 
lent boundary -layer mixing of such jets. 

The analysis has been carried out for the two classes of turbulence models: 

(1) methods based on differential -mean turbulence field without length scale equation 
(for example, the Donaldson-Gray model), and (2) methods based on differential -mean 
turbulence field with length scale equation (for example, the Spalding two -equation model). 
The methods of the first class are appropriate for near -equilibrium flows, whereas the 
methods of the second class are needed for analyzing the mixing of the nonequilibrium 
flows with differing scales. Predictions have been made by employing these models for 
a variety of reactive and nonreactive jets and comparison is made with selected experi- 
mental data. It is seen that the turbulence models (of the type of Spalding's two -equation 
model) with a length scale equation are adequate for most nonreactive flows of the type 
considered in this paper. Also, the Donaldson-Hilst invariant model needs some refine- 
ments of its "constants" to predict the axisymmetric jets analyzed here. 

For the reactive flows, where an allowance should be made for the concentration 
fluctuations in the finite -rate chemistry for an initially inhomogeneous mixture, an equa- 
tion similar to the concentration fluctuation equation of the Spalding three -equation model 
is suggested. This equation is much simpler to use than a similar proposal of Donaldson 
for an equilibrium field of turbulence. 

Because of the paucity of accurate experimental data for well-defined flow fields, 
some experimental data with simple chemistry such as 

N0 2 C1 + NO — * • NOC1 + N0 2 

are needed to test the equation proposed in this work to study the influence of the fluctua- 
tion (or mixedness) terms. 

Finally, the usefulness of the model proposed lies in handling the flows which have a 
chemistry classification (in relation to the molecular diffusion) in between the "fast" and 
"slow" regimes. These regimes have also been termed in the literature as "diffusion 
limited" and "reaction rate limited" regimes, respectively. 

Langley Research Center 

National Aeronautics and Space Administration 
Hampton, Va. 23665 
March 31, 1976 
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Figure 1. - Schematic of coaxial reactive jets. 



Figure 2.- Initial velocity profiles and characteristic shear 
widths for inner and outer regions of flow. 
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Figure 3.- Center-line velocity decay for krypton-air mixing 
with krypton injected centrally. 



Figure 4.- Center-line composition decay for krypton-air mixing 
with krypton injected centrally. 




o Experimental data from Forstall and Shapiro (ref. 15) 

Predictions from Donaldson-Gray model: 

: a = 1.0 



(a) Comparison with predictions from Donaldson-Gray model. 
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(b) Comparison with two-equation, Donaldson-Gray, 
and Donaldson -Hilst invariant models. 

Figure 5.- Center-line velocity decay for air -air mixing. 
34 
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Figure 6.- Axial distribution of velocity half-radius for air -air mixing. 
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Figure 7.- Center-line velocity decay for hydrogen-air mixing. 
U Q /U e = 3; Po/Pe = 0.08; \ = 0.24. 
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Figure 8.- Center-line composition decay for hydrogen-air mixing. 
U 0 /U e = 3; p Q /p e = 0.08; X = 0.24. 
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Figure 9.- Axial distribution of longitudinal turbulence intensity and 
center -line velocity decay for NO 2 CI-NO mixing with NO 2 CI 
injected centrally. 
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Donaldson-Gray model 

Two-equation model without cj^Cg in chemistry model 


Two -equation model with cj^Cg in chemistry model 

— Eddy -breakup model 



Figure 10.- Center -line composition decay for NO 2 CI-NO mixing 
with NO 2 CI injected centrally. 
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Figure 11.- Depletion rate of NO 2 CI at various axial locations. 
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